import math
import geojson

class MercatorProj:
    def __init__(self):
        self.iterative_times=10 #反解迭代次数
        self.iterative_value=0  #反解迭代初始值
        self.B0=0  #标准纬度（弧度制）
        self.L0=0  #原点经度（弧度制）
        self.A=6378137   #椭球体长半轴（米）,这里采用WGS84椭球体
        self.B=6356752.3142   #椭球体短半轴（米）

    def SetAB(self,a,b):
        """
        设定椭球体长短半轴
        """
        if a<=0 or b<=0:
            return
        self.A=a
        self.B=b


    def SetB0(self,b0):
        """
        设置标准纬度
        """
        if b0<-math.pi/2 or b0>math.pi/2:
            return
        self.B0=b0


    def SetL0(self,l0):
        """
        设置原点经度
        """
        if l0<-math.pi or l0>math.pi:
            return
        self.L0=l0


    def ToProj(self,L,B):
        """
        墨卡托投影正解，(B,L)->(X,Y)。因为geojson中节点存储形式为经度在前，纬度在后，
        故进行了调整
        """
        E=math.exp(1)

        if L<-math.pi or L>math.pi or B<-math.pi/2 or B>math.pi/2:
            return 1
        if self.A<=0 or self.B<=0:
            return 1
        
        f=(self.A-self.B)/self.A  #计算扁率

        dtemp=1-(self.B/self.A)*(self.B/self.A)
        if dtemp<0:
            return 1
        e1=math.sqrt(dtemp)  #计算第一偏心率

        dtemp=(self.A/self.B)*(self.A/self.B)-1
        if dtemp<0:
            return 1
        e2=math.sqrt(dtemp)  #计算第二偏心率

        #计算卯酉圈曲率半径N
        N=((self.A*self.A)/self.B)/math.sqrt(1+e2*e2*math.cos(self.B0)*math.cos(self.B0))
        K=N*math.cos(self.B0)
        #计算直角坐标系下的坐标
        Y=K*(L-self.L0)
        X=K*math.log(math.tan(math.pi/4+B/2)*pow((1-e1*math.sin(B))/(1+e1*math.sin(B)),e1/2))

        return X,Y
    

    def FromProj(self,X,Y):
        """
        墨卡托投影反解,(X,Y)->(B,L)
        """
        E=math.exp(1)

        if self.A<=0 or self.B<=0:
            return 1
        f=(self.A-self.B)/self.A #扁率

        dtemp=1-(self.B/self.A)*(self.B/self.A)
        if dtemp<0:
            return 1
        e1=math.sqrt(dtemp)  #第一偏心率

        dtemp=(self.A/self.B)*(self.A/self.B)-1
        if dtemp<0:
            return 1
        e2=math.sqrt(dtemp)  #第二偏心率

        N=((self.A*self.A)/self.B)/math.sqrt(1+e2*e2*math.cos(self.B0)*math.cos(self.B0))
        K=N*math.cos(self.B0)

        L=Y/K+self.L0

        #对B进行迭代
        B=self.iterative_value
        for i in range(self.iterative_times):
            B = math.pi / 2 - 2 * math.atan(pow(E, (-X / K)) * pow(E, (e1 / 2) * math.log((1 - e1 * math.sin(B)) / (1 + e1 * math.sin(B)))))

        return B,L

    @staticmethod
    def DegreeToRad(degree):
        """
        角度转弧度
        """
        return math.pi*(degree/180)
    
    @staticmethod
    def RadToDegree(rad1):
        """
        弧度转角度
        """
        return (rad1*180)/math.pi





def MercatorProjGeojson(input_file, output_file,is_forward):
    """
    从输入的GeoJSON文件中读取每个节点坐标，应用转换函数，
    然后将新坐标写入新的GeoJSON文件。

    输入参数：
    - input_file: 输入的GeoJSON文件路径
    - output_file: 输出的GeoJSON文件路径
    - is_forward: 定义正解或反解运算
    """
    with open(input_file, 'r',encoding='gb18030') as f:
        geojson_data = geojson.load(f)
        mercator_proj = MercatorProj()

    features = geojson_data['features']
    transformed_features = []

    for feature in features:
        geometry = feature['geometry']
        coordinates = geometry['coordinates']

        # 对每个节点坐标应用转换函数
        
        if is_forward==True:   #正解
            # transformed_coordinates = [MercatorProj.ToProj(MercatorProj.DegreeToRad(lat),MercatorProj.DegreeToRad(lon)) for lon, lat in coordinates]
            transformed_coordinates = [
        [mercator_proj.ToProj(MercatorProj.DegreeToRad(lon),
                              MercatorProj.DegreeToRad(lat)) for lon, lat in nested_coords]
        for nested_coords in coordinates
]
        else:  #反解
            transformed_coordinates=[]
           
            for nested_coords in coordinates:
                transformed_coords=[]
                for lon, lat in nested_coords:
                    new_lat=MercatorProj.RadToDegree(mercator_proj.FromProj(lon,lat)[0])  #转换后纬度，对应B
                    new_lon=MercatorProj.RadToDegree(mercator_proj.FromProj(lon,lat)[1])  #转换后经度，对应L
                    transformed_coords.append([new_lon, new_lat])
                transformed_coordinates.append(transformed_coords)

        # 更新新的坐标到GeoJSON feature
        transformed_feature = {'type': 'Feature',
                              'properties': feature['properties'],
                              'geometry': {'type': geometry['type'], 'coordinates': transformed_coordinates}}

        transformed_features.append(transformed_feature)

    transformed_geojson = {'type': 'FeatureCollection', 'features': transformed_features}
    with open(output_file, 'w',encoding='gb18030') as f:
        geojson.dump(transformed_geojson, f)
    print('done')